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We present an accelerated algorithm that samples correctly the thermodynamic ensemble in 
complex systems where the dynamics is controlled by activation barriers. The efficiency of the 
thermodynamically-weighted activation-relaxation technique (THWART) is many orders of magni- 
tude greater than standard molecular dynamics, even at room temperature and above, in systems 
as complex as proteins and amorphous silicon. 
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Throughout physics, chemistry and biology, a large 
proportion of atomistic processes take place on time 
scales many orders of magnitude longer than the typical 
phonon period. These processes are out of reach of single 
time-scale algorithms, such as molecular dynamics (MD), 
which can reach simulation times equivalent to the mi- 
crosecond at best. In view of this limitation, considerable 
effort has been devoted in the last few years to develop al- 
gorithms that allow spanning multiple time scales. These 
algorithms include the activation-relaxation technique 
(ART) 0, Q and similar techniques 0, 3 which sam- 
ple the configurational landscape by identifying transi- 
tion paths from minimum to minimum, as well as acceler- 
ated schemes based on molecular dynamics such as hyper- 
MD , temperature-assisted dynamics Q and others Q • 

ART and similar methods have been applied with suc- 
cess to study the topology of the energy landscape and 
activated mechanisms in a wide range of materials includ- 
ing amorphous and crystalline semiconductors H 0, 0] , 
glassy materials E2> atomic clusters 0, ^3 an d pro- 
teins El 111 El- ART defines events in the energy land- 
scape as a two-step process: (1) the system is first acti- 
vated from a local energy minimum to a nearby saddle- 
point and (2) then relaxed to a new minimum. Since 
the landscape constructed by ART consists only of local 
minima connected via saddle points, the entropic contri- 
butions to sampling are neglected so it is not possible to 
guarantee a proper statistical sampling, especially at el- 
evated temperatures where the harmonic approximation 
breaks down. 

Accelerated MD schemes, on the other hand, rely ei- 
ther on deforming the energy landscape [a, L| or on some 
projection from high-temperature simulations Q . In the 
first case, the energy basin surrounding a local minimum 
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is filled according to various rules. In the second case, 
multiple simulations are performed at high temperature 
and as soon as an event occurs, the time is rescaled with 
an appropriate factor. These methods work only for sim- 
ple systems, however, where it is possible to have a de- 
tailed a priori knowledge of the parameters defining the 
landscape. As soon as a wide range of barriers come 
into play, such as in proteins or disordered materials, it 
becomes very difficult to apply the appropriate transfor- 
mations in order to sample correctly and efficiently the 
phase space. 

In this Letter, we present an algorithm that provides 
a proper statistical sampling of the energy landscape at 
a wide range of temperatures irrespective of the com- 
plexity of the landscape. Combining molecular dynamics 
with ART, the thermodynamically-weighted activation- 
relaxation technique (THWART) samples the thermody- 
namically relevant parts of the phase space, hopping over 
barriers that can be many times higher than kuT. 

THWART generates motion in the energy landscape, 
a 3N-dimensional hypersurface, on which the height is 
given by the potential energy of a configuration as a 
function of its 3N atomic coordinates, where N is the 
number of atoms in the system. At low temperature, 
a configuration typically spends most of its time oscil- 
lating thermally around a local minimum, hopping over 
an energy barrier only when a thermal fluctuation trans- 
fers large amounts of energy onto a single mode. In this 
regime, it is possible to separate the energy landscape 
into two types of regions: basins and saddle regions. The 
basins are regions around local-energy minima where all 
Hessian eigenvalues, corresponding to the curvature of 
the landscape, are above a threshold value Ao- The sad- 
dle regions have at least one direction with a curvature 
(eigenvalue of the Hessian matrix) below this threshold 
value. These regions surround metastable points such a 
first or higher-order saddle points. A two-dimensional 
sketch of an energy landscape, paved according to these 
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FIG. 1: Sketch of a two-dimensional energy landscape. The 
black dots denote the locations of local-energy minima. These 
minima are part of basins, bounded by a line of constant 
lowest-curvature (solid line); the percolating region surround- 
ing the basins is called the saddle region. Basin-to-basin tra- 
jectories as generated by THWART are indicated by dashed 
lines. Constrained to ensure detailed balance, the trajectories 
come back to were they started if they fail to find a boundary. 



criteria, is shown in Fig. ^ 

At low temperature, the equilibrium properties of a 
system are determined by the basins, where the Boltz- 
mann weight exp(— E/ksT) is large and most of the 
sampling takes place. The partition function, Z = 
J dX exp(— E{X) /ksT), is thus well approximated by in- 
tegrating only over the basins. The basins, however, are 
disconnected regions in phase space and a proper sam- 
pling requires integrating over many of these basins. To 
ensure efficient sampling, it is therefore essential to ac- 
celerate the rate at which basins are visited. 

THWART achieves this acceleration in two steps. In 
the basin regions, configurations are evolved using stan- 
dard MD. We select the NVT ensemble and use the ve- 
locity Verlet algorithm for the integration. During the 
MD simulation, we monitor the lowest eigenvalue of the 
Hessian which defines the position of the boundary sep- 
arating the basin from the saddle region. As soon as this 
lowest eigenvalue reaches a given threshold Ao, the MD 
is suspended and the hopping phase starts. The atomic 
positions at the basin boundary is identified by xq and 
the velocities by vq. 

From xq, on the basin boundary, we construct a fully 
reversible path that leads to a new basin, crossing an en- 
ergy barrier. For this, we follow the eigendirection cor- 
responding to the lowest eigenvalue away from the basin 
until the eigenvalue crosses the threshold Ao from below. 
In order to ensure detailed balance, the trajectory into 
the saddle region is further constrained to move on a hy- 
perplane with constant potential energy. Specifically, the 
activated trajectory is generated by iterating the follow- 
ing equation: 



+ +c(fi;i + Fl;i+l) ■ (1) 
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Here, hi is the normalized eigenvector at x\ correspond- 



ing to the lowest (most negative) Hessian eigenvalue, F±. yi 
is the component of the force at Xi perpendicular to hi, 
At is a constant factor that determines the size of the 
increment, and c is a multiplicative constant, chosen to 
project the trajectory onto the hyperplane of constant 
potential energy. The orientation of ho is chosen initially 
so that it points towards the direction of more negative 
curvature, i.e., away from the initial basin; it is updated 
at each step by requiring that the inner product of the 
local eigenvector hi with that at the previous step, hi—i 
be always positive. 

The move along the eigendirection corresponding to 
the lowest eigenvalue changes the total configurational 
energy. The last term on the right-hand side of Eq. (JTJ 
corrects for this change and constrains the trajectory to 
the hyperplane with constant energy. Because the initial 
configuration is thermalized, its configurational energy 
is well above that of the local minimum, with roughly 
UbT/2 of additional potential energy per degree of free- 
dom. By transferring energy to a specific degree of free- 
dom from the heat bath formed by all the other 3A^ — 1 
degrees of freedom, the algorithm reproduces therefore 
the statistical mechanism responsible for crossing barri- 
ers. Simply moving along the force, to enforce this con- 
straint tends to bring the configuration back to its origi- 
nal position. Instead, we correct the energy by moving in 
the reduced space generated by the hyperplane perpen- 
dicular to this eigendirection, projecting the force onto 
this hyperplane, F±-i. 

Eq. (JIJ is iterated until the lowest eigenvalue passes the 
threshold (from below, this time) and the configuration 
reaches the boundary of a basin with position x p . At 
this point, the activation phase is stopped and the MD 
is resumed with the new positions and velocities v p = vq. 

The path generated from xq to x p is fully reversible: 
a configuration in basin p reaching x p would trigger the 
activation, bringing it to the other end of this path, in 
xq. Reversibility is ensured by the symmetric criterion 
for entering and leaving the saddle region as well as by 
keeping the path on a hyperplane of constant energy. Re- 
versibility is not sufficient for detailed balance, however. 
In continuous space, it is also necessary to verify that 
an infinitesimal volume surrounding xq remains constant 
as the configuration moves along the path leading to x p . 
This transformation, given by the Jacobian, determines 
the change in phase space, or entropy, between these two 
points. We have verified numerically that the Jacobian 
of transformation is unity along the whole path from one 
basin boundary to the other. In combination with re- 
versibility, the preservation of phase space ensures that 
detailed balance is respected. 

We also note that the path does not always lead to 
a new basin. In our simulations, it is not rare to see 
the trajectory form a circular path, coming back exactly 
at the initial point, xq, on the boundary after a long 
excursion in the saddle region. If the path does not close 
on itself, it generally connects to a different basin. 

Having established the validity of the algorithm, we 



apply THWART to two non-trivial systems: amorphous 
silicon and a small peptide. Simulations on amor- 
phous silicon are performed on models with two different 
sizes and initial configurations. The 1000-atom model 
was generated with ART nouveau Q using a modified 
Stillinger- Weber (mSW) potential fitted to the amor- 
phous phase 0| . This configuration is well-relaxed and 
is described in Ref. [l^. The 500-atom model was pro- 
duced with a bond-switching algorithm 0, 0] and sim- 
ply relaxed with the mSW potential, showing higher 
strain than the 1000-atom model. 

Both models were then evolved with MD in the stan- 
dard NVT ensemble at a temperature of 600 K and 800 
K, i.e., well below the melting transition temperature of 
^2000 K for this potential, but near the crystallization 
temperature from the amorphous phase. The MD simu- 
lations are integrated with time steps of 1 fs and run over 
a simulated time of 10 ns, except for the 500-atom con- 
figuration at 800 K, which is run for 20 ns. Figs. [21 and [3| 
show the total squared displacement, (r 2 ), between the 
initial configuration (which is energy-minimized) and the 
quenched configuration at time t, measured in the num- 
ber of force evaluations (with 1 force evaluation per time 
step), for these four distinct runs. 

For the well-relaxed 1000-atom model at 600 K, the fig- 
ure shows that after a few initial local rearrangements, 
the MD runs remain trapped around a few nearby min- 
ima for more than 9 ns. We find a similar situation for 
the 800 K simulation, with the configuration diffusing 
slightly further but still failing to explore more than a 
few nearby basins during the 10 ns simulation. 

The situation is quite different with THWART: In the 
basin, we follow the same MD procedure as above, while 
computing the lowest eigenvalue every 50 steps using a 
20-level Lanczos scheme. The MD run is continued until 
this lowest eigenvalue falls below the threshold value Ao- 
To avoid moving back and forth along the same path, the 
MD procedure is required to take at least 400 steps. From 
this point, we apply Eq. Q until the lowest eigenvalue 
increases above the same threshold. The total squared 
displacement measured as a function of the number of 
force evaluations for THWART is also plotted in Figs. [3 
and |3J 10 million force evaluations correspond roughly to 
5000 events. As can be seen, THWART explores the en- 
ergy landscape many orders of magnitude faster than MD 
at the same temperature. THWART is slightly slower at 
short times, because it performs one event at a time, 
while MD can activate multiples events in parallel across 
the model. However, it rapidly surpasses MD. In partic- 
ular, THWART does not seem to become trapped either 
at 600 K or 800 K. It is to be expected, however, that at 
high temperatures, when the number of events occurring 
in parallel starts proliferating, MD will become more effi- 
cient than THWART. From our simulations, this should 
happen close to the melting point. 

The efficiency of THWART depends on the value of 
the threshold, the sole parameter in THWART. Its value 
impacts the ratio of basin to saddle regions on the energy 
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FIG. 2: Total squared displacement as a function of the 
number of force evaluations, at temperatures 600 K (upper 
box) and 800 K (lower box), obtained with MD (circles) and 
THWART simulations (lines). For the THWART simulations 
at 600 K, we use A = -5A 2 (solid line) and -15 A 2 (dotted 
line). At 800 K, we use A = -5A 2 (solid line) and -10 A 2 
(dotted line). 



landscape; a very low value of Ao reduces THWART to 
MD, while a high value makes THWART more similar 
to ART. The impact of the threshold is shown in Fig. [21 
which plots the total squared displacement as a function 
of force evaluation for two values of Ao at 600 K and 800 
K. At both temperatures, THWART diffuses faster with 
a higher threshold, which reduces the basin size but also 
delivers a higher fraction of open paths. 

Even though the 500-atom model is more strained than 
the 1000-atom one, the configuration gets trapped very 
rapidly at 600 K with MD, within less than 1 ns. At 800 
K, however, the configuration is able to overcome many 
barriers and diffuse significantly, becoming trapped only 
after about 10 ns of simulation (see Fig. |3J). In spite of 
this considerable collective relaxation, THWART over- 
comes MD at 800 K, demonstrating the general efficiency 
of the method. 

Results on a small 10-residue peptide show that 
THWART is also more efficient than MD for this system. 
We use an artificial peptide of sequence AAAAAGAAAA 
with interactions described by CHARMM19 [H and an 
ASP solvation term [2lJ as implemented by Ponder in his 
program Tinker [22| . The peptide is first relaxed near its 
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FIG. 3: Squared displacement as a function of the number of 
force evaluations, at a temperature of 800 K, obtained with 
MD (circles) and THWART simulations (lines) with Ao = 
— 20A 2 . The inset shows the same at 600 K, and with Ao = 
-5A 2 . 




FIG. 4: Superimposed energy-minimized backbone configu- 
rations of the peptide at different times, in MD (left) and 
THWART (right) simulations at a temperature of 600K. 
About 100 configurations are superimposed for each method. 



energy minimum using THWART. This sequence pos- 
sesses a large number of metastable minima surrounding 
this lowest-energy state. 

We first perform a 300-ps constant-temperature MD 
simulation using Tinker. To characterize the efficiency of 
sampling, the energy-minimized conformations after ev- 
ery 2 ps are graphically superimposed in the left panel 
of Fig. ^ Clearly, the conformation is trapped in a few 
nearby local minima. Next, a THWART simulation is 
performed of the same system over 300 000 force evalu- 
ations, or about 150 events. The energy-minimized con- 
formations are superimposed in the right panel of Fig.0] 
THWART clearly samples the phase space around the 
initial state much more efficiently than MD. 

In summary, the thcrmodynamically-weighted 
activation-relaxation technique is an accelerated algo- 
rithm for sampling efficiently the relevant parts of the 
phase space in systems where the dynamics is controlled 
by activation barriers. THWART samples correctly the 
thermodynamic ensemble, with an efficiency many orders 
of magnitude greater than standard molecular dynamics, 
even at room temperature and above. We have shown 
that THWART is efficient in various complex systems 
such as proteins and amorphous silicon. THWART 
should be particularly useful to study materials such as 
glasses and proteins in water, where other accelerated 
techniques fail. 
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